#include <math.h>
#include <stdlib.h>
#include <string.h>
#include <stdio.h>
#include <stdbool.h>
#include "Math.h"

/* Math library created by Daniel White
	Various math functions that could be useful for game developement (and then some more)
*/

/* Globals */
int *Math_PrimeCache=NULL;
int Math_PrimeCacheSize=0;

/* Functions */

int Math_GetFirstNPrimes(int Amount, int *Array)
{
	if (Amount<=Math_PrimeCacheSize)
	{
		memcpy(Array, Math_PrimeCache, Amount);
		return 0;
	}
	
	int UpdateCache, *TempPointer=realloc(Math_PrimeCache, (Amount+1)*sizeof(Math_PrimeCache));
	
	if (TempPointer==NULL)
		UpdateCache=0;
	else
	{
		UpdateCache=1;
		Math_PrimeCacheSize=Amount;
		Math_PrimeCache=TempPointer;
	}
	
	/* Work out an upper bound for the Amount'th prime number */
	double A=2, B=4;
	while(A!=B)
	{
		if (((int)(B/log(B)))<Amount)
		{
			A=B;
			B*=1.15;
		}
		if (((int)(B/log(B)))>Amount)
			B*=0.85;
		if (((int)(B/log(B)))==Amount)
			A=B;
	}
	long Max=(long)A;
	
	/* Now sieve */
	long X=3, X2, Y, AmountFound=0, SqrtMax=sqrt(Max);
	
	bool *Numbers=calloc(Max, sizeof(bool));
	if (Numbers==NULL)
		return 1;
	
	Array[0]=2;
	if (UpdateCache)
		Math_PrimeCache[0]=2;
	
	while(X<SqrtMax)
	{
		if (Numbers[X]!=1)
		{
			++AmountFound;
			if (UpdateCache)
				Math_PrimeCache[AmountFound]=X;
			Array[AmountFound]=X;
			
			X2=X<<1;
			for(Y=X;Y<Max;Y+=X2)
				Numbers[Y]=1;
		}
		X+=2;
	}
	
	while(AmountFound<Amount)
	{
		if (Numbers[X]!=1)
		{
			++AmountFound;
			if (UpdateCache)
				Math_PrimeCache[AmountFound]=X;
			Array[AmountFound]=X;
		}
		X+=2;
	}
	
	free(Numbers);
	
	return 0;
}

int Math_GetFirstNPrimesLowMemory(int Amount, int *Array)
{	
	if (Amount<=Math_PrimeCacheSize)
	{
		memcpy(Array, Math_PrimeCache, Amount);
		return 0;
	}
	
	Array[0]=2;
	Array[1]=3;
	int Number=5, Found=2, Divisor, SqrtNumber, IsPrime;
	for(;Found<Amount;)
	{
		SqrtNumber=sqrt(Number);
		IsPrime=1;
		for(Divisor=3;Divisor<=SqrtNumber;Divisor+=2)
		{
			if (Number%Divisor==0)
			{
				IsPrime=0;
				Divisor=SqrtNumber;
			}
		}
		if (IsPrime)
		{
			Array[Found]=Number;
			++Found;
		}
		
		Number+=2;
		if (Found==Amount)
			return 0;
			
		SqrtNumber=sqrt(Number);
		IsPrime=1;
		for(Divisor=3;Divisor<=SqrtNumber;Divisor+=2)
		{
			if (Number%Divisor==0)
			{
				IsPrime=0;
				Divisor=SqrtNumber;
			}
		}
		if (IsPrime)
		{
			Array[Found]=Number;
			++Found;
		}
			
		Number+=4;
	}
	
	return 0;
}

int Math_IsPrime(int Number)
{
	int Amount=100;
	int TempArray[Amount];
	if (Math_PrimeCacheSize<Amount)
		Math_GetFirstNPrimes(Amount, TempArray);
	
	/* Divide by primes */
	int PrimeIndex;
	for(PrimeIndex=0;PrimeIndex<Math_PrimeCacheSize;++PrimeIndex)
	{
		if (Number==Math_PrimeCache[PrimeIndex])
			return 1;
		else
		{
			while (Number%Math_PrimeCache[PrimeIndex]==0)
				Number/=Math_PrimeCache[PrimeIndex];
		}
	}
	
	if (Number==1)
		return 0;
	
	/* Now try trial division */
	int SqrtNumber=sqrt(Number), Divisor;
	for(Divisor=3;Divisor<=SqrtNumber;Divisor+=2)
	{
		if (Number%Divisor==0)
			return 0;
	}
	
	return 1;
}

int Math_GCD(int Number0, int Number1)
{
	int Temp;
	while(Number1!=0)
	{
		Temp=Number1;
		Number1=Number0%Number1;
		Number0=Temp;
	}
	return Number0;
}

inline int Math_LCM(int Number0, int Number1)
{
	return abs(Number0*Number1)/Math_GCD(Number0, Number1);
}

int Math_DigitalSum(int Number)
{
	int Sum=0;
	
	while(Number>0)
	{
		Sum+=Number%10;
		Number/=10;
	}
	
	return Sum;
}

int Math_DigitalRoot(int Number)
{
	if (Number%9==0)
		return 9;
	else
		return Number%9;
}

inline float Math_AngleY(float Angle)
{
	if(Angle < 90.0)
		return -cos((float)Math_DegToRad(Angle));
	if(Angle < 180.0)
		return sin((float)Math_DegToRad(Angle-90.0));
	if(Angle < 270.0)
		return cos((float)Math_DegToRad(Angle-180.0));
	return -sin((float)Math_DegToRad(Angle-270.0));
}

inline float Math_AngleX(float Angle)
{
	if(Angle < 90.0)
		return sin((float)Math_DegToRad(Angle));
	if(Angle < 180.0)
		return cos((float)Math_DegToRad(Angle-90.0));
	if(Angle < 270.0)
		return -sin((float)Math_DegToRad(Angle-180.0));
	return -cos((float)Math_DegToRad(Angle-270.0));
}

inline float Math_XYAngle(float X1, float Y1, float X2, float Y2)
{
	return Math_RadToDeg(atan2(Y2-Y1, X2-X1))+90;
}

float Math_NormaliseAngle(float Angle)
{
	while(Angle >= 360.0)
		Angle -= 360.0;
	while(Angle < 0.0)
		Angle += 360.0;
	return Angle;
}

int Math_PrimeFactors(int Number, int **Factors)
{
	int TempPrimes[10];
	Math_GetFirstNPrimes(10, TempPrimes);
	
	int FactorCount=0, PrimeIndex, *TempPointer;
	
	for(PrimeIndex=0;PrimeIndex<Math_PrimeCacheSize;++PrimeIndex)
	{
		while (Number%Math_PrimeCache[PrimeIndex]==0)
		{
			Number/=Math_PrimeCache[PrimeIndex];
			++FactorCount;
			
			TempPointer=realloc((*Factors), FactorCount*sizeof(*Factors));
			if (TempPointer==NULL)
				return -1;
			(*Factors)=TempPointer;
			(*Factors)[FactorCount-1]=Math_PrimeCache[PrimeIndex];
		}
	}
	
	int Divisor;
	for(Divisor=3;Divisor<=Number;Divisor+=2)
	{
		while (Number%Divisor==0)
		{
			Number/=Divisor;
			++FactorCount;
			
			TempPointer=realloc((*Factors), FactorCount*sizeof(*Factors));
			if (TempPointer==NULL)
				return -1;
			(*Factors)=TempPointer;
			(*Factors)[FactorCount-1]=Divisor;			
		}
		
		if (Number==1)
			return FactorCount;
	}
	
	return FactorCount;
}

int Math_PrimeFactorCount(int Number)
{
	int *PrimeFactors=NULL;
	int PrimeAmount=Math_PrimeFactors(Number, &PrimeFactors);
	free(PrimeFactors);	
	return PrimeAmount;
}

int Math_FactorCount(int Number)
{
	int *PrimeFactors=NULL;
	int PrimeAmount=Math_PrimeFactors(Number, &PrimeFactors);
	if (PrimeAmount==-1)
		return -1;
	
	int Amount=1, Last=0, Index, Count=1;
	for(Index=0;Index<PrimeAmount;++Index)
	{
		if (PrimeFactors[Index]!=Last)
		{
			Amount*=Count;
			
			Last=PrimeFactors[Index];
			Count=2;
		}
		else
			++Count;
	}
	Amount*=Count;
	
	free(PrimeFactors);	
	return Amount;
}

void Math_ClearCache()
{
	/** if (Math_PrimeCache!=NULL) */
		free(Math_PrimeCache);
	Math_PrimeCache=NULL;
	Math_PrimeCacheSize=0;
}
